/******************************************************************************/
/*                                                                            */
/*  Copyright (c) 2007-2009: ngi-central.org                                  */
/*  All rights reserved.                                                      */
/*                                                                            */
/******************************************************************************/


/// \file
/// This file implements data structures for moments calculations.


#ifndef __sgl_central_moments_h__
#define __sgl_central_moments_h__

#include <sgl/global_settings.h>
#include <sgl/conditions.h>
#include <sgl/normalized_moments.h>
#include <sgl/real_concept.h>
#include <boost/concept_check.hpp>

namespace sgl {


/******************************************************************************/
/*                                                                            */
/* central_moments                                                            */
/*                                                                            */
/******************************************************************************/

template <
	/// The type. 
	//! 
	//! This type defaults to SGL_FPTYPE (which in turn defaults to double), but 
	//! any other numeric type is also suitable.
	typename T=SGL_FPTYPE>
/// Class to hold central moments up to the third order.
//!
//! Central moments are invariant with respect to translation, but not 
//! invariant with respect to scaling and rotation.
//!
//! The following operations are implemented:
//! operator== : comparison for equality, directly implemented.
//! operator != : comparison for inequality, implemented through 
//! boost::equality_comaparable.
class central_moments
	: boost::equality_comparable<sgl::central_moments<T>, 
		sgl::central_moments<T>
	>
{
	// Only real types qualify for type T. If this line triggers an error,
	// type T is not a real type.
	BOOST_CONCEPT_ASSERT((sgl::real_concept<T>));

public:
	/// The value type.
	typedef T value_type;

	/// Default constructor. 
	//!
	//! Zero initializes the moments.
	central_moments()
		: mu00_(T())
		, mu20_(T())
		, mu11_(T())
		, mu02_(T())
		, mu30_(T())
		, mu21_(T())
		, mu12_(T())
		, mu03_(T())
	{ }

	/// Standard constructor.
	//!
	//! Initializes the moments.
	central_moments(
		/// Zeroth order moment.
		T mu00,
		/// Second order horizontal moment.
		T mu20, 
		/// Second order mixed moment.
		T mu11,
		/// Second order vertical moment.
		T mu02,
		/// Third order horizontal moment.
		T mu30, 
		/// Third order horizontal mixed moment.
		T mu21,
		/// Third order vertical mixed moment.
		T mu12,
		/// Third order vertical moment.
		T mu03)
		: mu00_(mu00)
		, mu20_(mu20)
		, mu11_(mu11)
		, mu02_(mu02)
		, mu30_(mu30)
		, mu21_(mu21)
		, mu12_(mu12)
		, mu03_(mu03)
	{ }

	/// Constructs central moments from normalized moments.
	//!
	//! This calculates central moments from normalized moments.
	central_moments(
		/// The right hand side.
		sgl::normalized_moments<T> const& r)
		: mu00_(r.calculate_mu00())
		, mu20_(r.calculate_mu20())
		, mu11_(r.calculate_mu11())
		, mu02_(r.calculate_mu02())
		, mu30_(r.calculate_mu30())
		, mu21_(r.calculate_mu21())
		, mu12_(r.calculate_mu12())
		, mu03_(r.calculate_mu03())
	{ }

	/// Constructs a copy.
	central_moments(
		/// The right hand side.
		sgl::central_moments<T> const& r)
		: mu00_(r.mu00_)
		, mu20_(r.mu20_)
		, mu11_(r.mu11_)
		, mu02_(r.mu02_)
		, mu30_(r.mu30_)
		, mu21_(r.mu21_)
		, mu12_(r.mu12_)
		, mu03_(r.mu03_)
	{ }

	/// Copy assignment
	sgl::central_moments<T> & operator =(
		/// The right hand side.
		sgl::central_moments<T> rhs)
	{
		swap(rhs);
		return *this;
	}

	/// Exchanges two objects.
	void swap(
		/// The other object.
		sgl::central_moments<T> & rhs)
	{
		std::swap(mu00_, rhs.mu00_);
		std::swap(mu20_, rhs.mu20_);
		std::swap(mu11_, rhs.mu11_);
		std::swap(mu02_, rhs.mu02_);
		std::swap(mu30_, rhs.mu30_);
		std::swap(mu21_, rhs.mu21_);
		std::swap(mu12_, rhs.mu12_);
		std::swap(mu03_, rhs.mu03_);
	}

	/// Compare objects for equality.
	//!
	//! This instance is equal to another instance, if its moments are all equal 
	//! to the moments of the other instance.
	//!
	//! Through the use of boost::equality_comparable operator != is also
	//! defined and implemented by calling this operator.
	//!
	//! \return true if equal, false otherwise.
    bool operator ==(
		/// Right hand side.
		sgl::central_moments<T> const& rhs) const
    {
        return 
			mu00_ == rhs.mu00_ && 
			mu20_ == rhs.mu20_ && 
			mu11_ == rhs.mu11_ &&
			mu02_ == rhs.mu02_ && 
			mu30_ == rhs.mu30_ && 
			mu21_ == rhs.mu21_ &&
			mu12_ == rhs.mu12_ &&
			mu03_ == rhs.mu03_;
    }

	/// Access the zeroth order central moment.
	//!
	//! /return Returns a reference to the zeroth order central moment.
	T & mu00()
	{
		return mu00_;
	}

	/// Access the zeroth order central moment.
	//!
	//! /return Returns a const reference to the zeroth order central moment.
	T const & mu00() const
	{
		return mu00_;
	}

	/// Access the horizontal second order central moment.
	//!
	//! /return Returns a reference to the horizontal second order central 
	//! moment.
	T & mu20()
	{
		return mu20_;
	}

	/// Access the horizontal second order central moment.
	//!
	//! /return Returns a const reference to the horizontal second order central 
	//! moment.
	T const & mu20() const
	{
		return mu20_;
	}

	/// Access the mixed second order central moment.
	//!
	//! /return Returns a reference to the mixed second order central moment.
	T & mu11()
	{
		return mu11_;
	}

	/// Access the mixed second order central moment.
	//!
	//! /return Returns a const reference to the mixed second order central 
	//! moment.
	T const & mu11() const
	{
		return mu11_;
	}

	/// Access the vertical second order central moment.
	//!
	//! /return Returns a reference to the vertical second order central moment.
	T & mu02()
	{
		return mu02_;
	}

	/// Access the vertical second order central moment.
	//!
	//! /return Returns a const reference to the vertical second order central 
	//! moment.
	T const & mu02() const
	{
		return mu02_;
	}

	/// Access the horizontal third order central moment.
	//!
	//! /return Returns a reference to the horizontal third order central 
	//! moment.
	T & mu30()
	{
		return mu30_;
	}

	/// Access the horizontal third order central moment.
	//!
	//! /return Returns a const reference to the horizontal third order central 
	//! moment.
	T const & mu30() const
	{
		return mu30_;
	}

	/// Access the first (horizontal) mixed third order central moment.
	//!
	//! /return Returns a reference to the first (horizontal) mixed third order 
	//! central moment.
	T & mu21()
	{
		return mu21_;
	}

	/// Access the first (horizontal) mixed third order central moment.
	//!
	//! /return Returns a const reference to the first (horizontal) mixed third 
	//! order central moment.
	T const & mu21() const
	{
		return mu21_;
	}

	/// Access the second (vertical) mixed third order central moment.
	//!
	//! /return Returns a reference to the second (vertical) mixed third order 
	//! central moment.
	T & mu12()
	{
		return mu12_;
	}

	/// Access the second (vertical) mixed third order central moment.
	//!
	//! /return Returns a const reference to the second (vertical) mixed third 
	//! order central moment.
	T const & mu12() const
	{
		return mu12_;
	}

	/// Access the vertical third order central moment.
	//!
	//! /return Returns a reference to the vertical third order central moment.
	T & mu03()
	{
		return mu03_;
	}

	/// Access the vertical third order central moment.
	//!
	//! /return Returns a const reference to the vertical third order central 
	//! moment.
	T const & mu03() const
	{
		return mu03_;
	}

	/// Calculate eta20 from central moments.
	//!
	//! eta20 is the second order horizontal scale invariant moment.
	//!
	//! /return Second order horizontal scale invariant moment.
	T calculate_eta20() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		return mu20() / (mu00()*mu00());
	}

	/// Calculate eta11 from central moments.
	//!
	//! eta11 is the second order mixed scale invariant moment.
	//!
	//! /return Second order mixed scale invariant moment.
	T calculate_eta11() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		return mu11() / (mu00()*mu00());
	}

	/// Calculate eta02 from central moments.
	//!
	//! eta02 is the second order vertical scale invariant moment.
	//!
	//! /return Second order vertical scale invariant moment.
	T calculate_eta02() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		return mu02() / (mu00()*mu00());
	}

	/// Calculate eta30 from central moments.
	//!
	//! eta30 is the third order horizontal scale invariant moment.
	//!
	//! /return Third order horizontal scale invariant moment.
	T calculate_eta30() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		return mu30() / pow(mu00(), 5./2.);
	}

	/// Calculate eta21 from central moments.
	//!
	//! eta21 is the third order horizontal mixed scale invariant moment.
	//!
	//! /return Third order horizontal mixed scale invariant moment.
	T calculate_eta21() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		return mu21() / pow(mu00(), 5./2.);
	}

	/// Calculate eta12 from central moments.
	//!
	//! eta12 is the third order vertical mixed scale invariant moment.
	//!
	//! /return Third order vertical mixed scale invariant moment.
	T calculate_eta12() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		return mu12() / pow(mu00(), 5./2.);
	}

	/// Calculate eta03 from central moments.
	//!
	//! eta03 is the third order vertical scale invariant moment.
	//!
	//! /return Third order vertical scale invariant moment.
	T calculate_eta03() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		return mu03() / pow(mu00(), 5./2.);
	}

	/// Calculate the first Flusser moment from central moments.
	//!
	//! fm1 is the first scale and rotation invariant Flusser moment.
	//!
	//! /return First scale and rotation invariant Flusser moment.
	T calculate_fm1() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		T a = mu20()*mu02();
		T b = mu11()*mu11();
		T c = mu00()*mu00();
		return (a - b) / (c*c);
	}

	/// Calculate the second Flusser moment from central moments.
	//!
	//! phi2 is the second scale and rotation invariant Flusser moment.
	//!
	//! /return Second scale and rotation invariant Flusser moment.
	T calculate_fm2() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		T a = mu30()*mu30();
		T b = mu03()*mu03();
		T c = mu30()*mu21()*mu12()*mu03();
		T d = mu30()*mu12()*mu12()*mu12();
		T e = mu03()*mu21()*mu21()*mu21();
		T f = mu21()*mu21();
		T g = mu12()*mu12();
		T h = mu00()*mu00()*mu00();
		return (a*b - T(6)*c + T(4)*(d + e) - T(3)*f*g)/(h*h*h*mu00());
	}

	/// Calculate the third Flusser moment from central moments.
	//!
	//! phi3 is the third scale and rotation invariant Flusser moment.
	//!
	//! /return Third scale and rotation invariant Flusser moment.
	T calculate_fm3() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		T a = mu21()*mu03();
		T b = mu12()*mu12();
		T c = mu30()*mu03();
		T d = mu21()*mu12();
		T e = mu30()*mu12();
		T f = mu21()*mu21();
		T g = mu20()*(a - b);
		T h = mu11()*(c - d);
		T i = mu02()*(e - f);
		T j = mu00()*mu00()*mu00();
		return (g - h + i)/(j*j*mu00());
	}

	/// Calculate the fourth Flusser moment from central moments.
	//!
	//! phi4 is the fourth scale and rotation invariant Flusser moment.
	//!
	//! /return Fourth scale and rotation invariant Flusser moment.
	T calculate_fm4() const
	{
		SGL_PRECONDITION(mu00() != T(0), "division by zero.");

		T a = mu20()*mu20()*mu20()*mu03()*mu03();
		T b = mu20()*mu20()*mu11()*mu12()*mu03();
		T c = mu20()*mu20()*mu02()*mu21()*mu03();
		T d = mu20()*mu20()*mu02()*mu12()*mu12();
		T e = mu20()*mu11()*mu11()*mu21()*mu03();
		T f = mu20()*mu11()*mu02()*mu30()*mu03();
		T g = mu20()*mu11()*mu02()*mu21()*mu12();
		T h = mu11()*mu11()*mu11()*mu30()*mu03();
		T i = mu20()*mu02()*mu02()*mu30()*mu12();
		T j = mu20()*mu02()*mu02()*mu21()*mu21();
		T k = mu11()*mu11()*mu02()*mu30()*mu12();
		T l = mu11()*mu02()*mu02()*mu30()*mu21();
		T m = mu02()*mu02()*mu02()*mu30()*mu30();
		T n = mu00()*mu00()*mu00()*mu00()*mu00();
		return (a + m + - T(8)*h + T(3)*(T(2)*(-b - c + f - i - l) + 
			T(3)*(d + j) + T(4)*(e + k) - T(6)*g))/(n*n*mu00());
	}

private:
	/// The zeroth order central moment.
	T mu00_;

	/// The horizontal second order central moment.
	T mu20_;

	/// The mixed second order central moment.
	T mu11_;

	/// The vertical second order central moment.
	T mu02_;

	/// The horizontal third order central moment.
	T mu30_;

	/// The first (horizontal) mixed third order central moment.
	T mu21_;

	/// The second (vertical) mixed third order central moment.
	T mu12_;

	/// The vertical third order central moment.
	T mu03_;
};


} // namespace sgl

#endif // __sgl_central_moments_h__
